function [aa,maxB]=magnet4(v0,L,C,R,U0,l,lcoil)
%v0 is the initial velocity;
%L is inductor
%C is capacity
%R is resistence
%R^2<4LC must be satisfied
%U0 is the initial volumn
%l is the length of the thing
%lcoil is the distance between two coils
%neglect all the wave of currency after t=T/2.
%a typical input is magnet4(100,0.000001,0.01,0.0001,1000000,1,0.1)
%another typical input is magnet4(100,0.0001,0.1,0.005,5000000,1,3)

miu0=4*pi*1e-7;
r=sqrt(1/pi);%the radius of the thing
A=pi*r^2;%A is the area of the thing
rcoil=1;%rcoil is the radius of the coil
w=sqrt(4*L*C-R^2)/(2*L*C);
TT=2*pi/w;
Bm=1.4;%the magnet of the thing
m=Bm/miu0*A*sqrt(l^2+(2*r)^2);%m is the ciju..
M=10000;%the mass of the thing

    function Itans=It(t)
        Itans=2*C*exp(-R*t/(2*L*C))*U0.*sinh(sqrt(-4*L*C+R^2)*t/(2*L*C))/sqrt(-4*L*C+R^2);
    end

    function Bxtans=Bxt(x,t)
        Bxtans=miu0*rcoil*rcoil*It(t)/2/(rcoil*rcoil+x*x)^(3/2);
    end

%when the top of the thing reaches the position of x=0, the currency needs
%to be maximum. That means that when the centre of the thing is at the
%point of x=-l/2-T/4*v0, we turn the capacity on, and we call that moment t=0.

    function dy=f(t,y)
        dy=zeros(2,1);%y(1)=h;y(2)=v;y(3)=a;
        dy(1)=y(2);
        dy(2)=m/M/l*(Bxt(y(1)+l/2,t)-Bxt(y(1)-l/2,t));
    end

h0=-l/2-TT/4*v0;
options = odeset('RelTol',1e-9,'AbsTol',1e-9);    
[T,Y]=ode15s(@f,[0 0.1],[h0 v0],options);
lines=size(Y(:,1),1);
n=lines;
for ii=1:lines
    if Y(ii,1)>(h0+lcoil)
        n=ii;
        break;
    end;
end;

h=zeros(n,1);
v=zeros(n,1);
t=zeros(n,1);
a=zeros(n,1);
inta=0;
for ii=1:n
    h(ii)=Y(ii,1);
    v(ii)=Y(ii,2);
    t(ii)=T(ii);
    a(ii)=m/M/l*(Bxt(h(ii)+l/2,t(ii))-Bxt(h(ii)-l/2,t(ii)));
    if ii>1
        inta=inta+a(ii)*(h(ii)-h(ii-1));
    end;
end;

maxB=0;
for iii=1:n
    if (Bxt(0,t(iii))>maxB) 
        maxB=Bxt(0,t(iii));
    end
end

tt=t(n);
aa=inta/(h(n)-h(1));
 figure(1);
 plot(h,a); 
 title('a~h graph');
 xlabel('h/m');
 ylabel('a/(m/s^2)');
 figure(2);
 plot(t,a);
 title('a~t graph');
 xlabel('t/s');
 ylabel('a/(m/s^2)');
 figure(3);
 plot(t,Bxt(0,t));
 title('B~t graph');
 xlabel('t/s');
 ylabel('B/T');

end